Refactoring

Following #3411 it is clear that the grids require some refactor to allow ZStar.
This PR implements another proposal for ZStar which hinges on changing the grids to have just one z field.
The case of a static z (the only allowed in Oceananigans at the moment) is covered with a StaticVerticalCoordinate type that contains all the specifications that were in the grid before:

This refactor should change nothing in the user interface nor the internals of the code, since the whole API should depend on znode and Δz, but just reorganizes the vertical coordinate in a compact type that can be modified.

ZStar Implementation

With this refactor, a ZStar coordinate is implemented through a ZStarVerticalCoordinate type

e₃ is the vertical scaling of the grid spacing defined as

$$\frac{\partial z}{\partial r} = \frac{H + \eta}{H}$$

with $z$ the actual spatial coordinate (moving with the free surface) and $r$ is the reference vertical coordinate also called $z^\star$ (equivalent to a static coordinate, or $z$ when $\eta = 0$).
ηⁿ is the value of the free surface at the current time step, needed to calculate znodes

$$z = r \cdot e_3 + \eta \ .$$

Finally $H$ is the output of static_column_depth.
As a consequence of the definitions, all $z$ properties (znodes, zspacings, Δz, etc...) are replaced with $r$ (rnodes, rspacings,Δr, etc...) and the $z$ properties for a zstar vertical coordinate are calculated appropriately in the Operators module (https://github.com/CliMA/Oceananigans.jl/blob/ss/new-zstar/src/Operators/variable_grid_operators.jl)

An example of the proposed user interface is shown in the validation/z_star_coordinate examples

Changes in the model timestepping

Another requirement is changing the internals of the HydrostaticFreeSurfaceModel to solve the correct equations, in particular, the notable changes are:

(1) Updating the grid at each timestep in the update_state! step

(2) including the gradient of the grid in the momentum equations, calculated as:

(3) Changing the computation of the vertical velocity to include the movement of the grid

(4) Advancing $e_3\theta$ instead of $\theta$ in the tracer equations and subsequently unscale it back after the grid scaling factor at the new time step $e_3^{n+1}$ is known

(5) Adding a non-linear free surface by changing the static_column_depth in the split explicit free surface implementation to a dynamic_column_depth defined as

Note: these changes are valid only for the QuasiAdamsBashort2 timestepper. Support for the SplitRungeKutta3 timestepper requires a bit more work and can be done when all the infrastructure is in place

OutputWriters?

The last piece of the puzzle (not implemented in this PR) would be to change the OutputWriters to include by default the variable grid properties in the timeseries field of the jld2 / netcdf writer.

Possibility for improvements

I got a bit carried away and completed a working implementation to make sure it could be done this way, but I am happy to change just about everything in here. I would like to know what people think about this implementation and what people suggest especially in terms of
(1) Implementation
(2) Variable naming

I can also split this PR in a couple of ones, probably one that includes the refactor to the grids (implementation of a StaticVerticalCoordinate) and one that implements a working version of a ZStarVerticalCoordinate.